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On the thermo-elastostatics of heterogeneous materials 
I. General integral equation 



Abstract We consider a linearly thermoelastic composite medium, which consists of a homogeneous matrix 
containing a statistically inhomogeneous random set of inclusions, when the concentration of the inclusions is 
a function of the coordinates (so-called Functionally Graded Materials) . The composite medium is subjected to 
essentially inhomogeneous loading by the fields of the stresses, temperature and body forces (e.g. for a centrifugal 
load). The general integral equations connecting the stress and strain fields in the point being considered and the 
surrounding points are obtained for the random fields of inclusions. The method is based on a centering procedure 
of subtraction from both sides of a known initial integral equation their statistical averages obtained without any 
auxiliary assumptions such as, e.g., effective field hypothesis implicitly exploited in the known centering methods. 
In so doing the size of a region including the inclusions acting on a separate one is finite, i.e. the locality principle 
takes place. 

Keywords: A. microstructures, B. inhomogeneous material, B. elastic material. 



1. Introduction 

The need for consideration of the actual microstructure of composite materials subjected to essentially 
inhomogeneous mechanical, body force and temperature loading in micromechanics problems is well 
known. Unfortunately, the starting assumptions made in the majority of studies, namely that the structure 
of the composite media as well as the random fields of stresses are statistically homogeneous and therefore 
are invariant with respect to the translation may be invalid. 

For example, due to some production technologies, the inclusion concentration may be a function of 
the coordinates (see e.g. [1-3]). The accumulation of damage also occurs locally in stress-concentration 
regions, for example, at the tip of a macroscopic crack (see e.g. [4]). Furthermore, in layered composite 
shells the location of the fibers is random within the periodic layers, and the micromechanics equations 
are equations with almost periodic coefficients. Finally, Functionally Graded Materials (FGMs) have been 
the subject of intense research efforts from the mid-1980s when this term was originated in Japan in the 
framework of a national project to develop heat-shielding structural materials for the future Japanese 
space program. FGM is a composite consisting of two or more phases which is fabricated with a spatial 
variation of its composition that may improve the structural response (see e.g. [5, 6])Q. Moreover, modern 
constructions from composite materials on frequent occasions are subjected to essentially inhomogeneous 
loading by fields of the stresses, temperature and body forces (e.g. for a centrifugal load). 

The final goals of micromechanical research of composites involved in a prediction of both the overall 
effective properties and statistical moments of stress-strain fields are based on the approximate solution 
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T A popular macroscopic approach is the modeling of FGMs as macroscopically elastic materials, in which the material 
properties are graded but continuous and are described by a local constitutive equation (see e.g. [7, 8]), Nevertheless, 
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detected in [9] (see also Chapter 12 in [10]). However, the problem of estimation of effective properties of such materials is 
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of exact initial integral equations connecting the random stress fields at the point being considered and 
the surrounding points. This infinite system of coupled integral equations is well-known for statistically 
homogeneous composite materials subjected to homogeneous boundary conditions (see e.g. [10-13]). The 
goal of this paper is to obtain a generalization of these equations for the case of statistically inhomogeneous 
structures of composite materials subjected to essentially inhomogeneous loading by fields of the stresses, 
temperature and body forces. The method is based on a centering procedure of subtraction from both sides 
of a known initial integral equation the statistical averages obtained without any auxiliary assumptions 
such as, e.g., effective field hypothesis implicitly exploited in the known centering methods. Working 
with statistical averages (rather than with volume averages) is convenient because statistical averaging 
commutes with differentiating and integrating that becomes fundamentally important for statistically 
inhomogeneous media. 

2. Preliminaries 

2.1 Basic equations 

Let a linear elastic body occupy an open simply connected bounded domain w C R d with a smooth 
boundary T and with an indicator function W and space dimensionality d (d = 2 and d = 3 for 2-D 
and 3-D problems, respectively). The domain w contains a homogeneous matrix and a statistically 
inhomogeneous set X = of inclusions vi with indicator functions V% and bounded by the closed 
smooth surfaces Ti (i = 1,2,...). It is assumed that the inclusions can be grouped into components 
(phases) (q = 1, 2, . . . , N) with identical mechanical and geometrical properties (such as the shape, 

size, orientation, and microstructure of inclusions). For the sake of definiteness, in the 2-D case we will 
consider a plane-strain problem. At first no restrictions arc imposed on the elastic symmetry of the phases 
or on the geometry of the inclusions, d The local strain tensor e is related to the displacements u via the 
linearized strain-displacement equation e = \ [V ® u + ( V <£> u) T ]. Here ® denotes tensor product, and 
(.) T denotes matrix transposition. The stress tensor, <x, satisfies the equilibrium equation:V • <x = — f, 
where the body force tensor f can be generated, e.g., by either gravitational loads or a centrifugal load. 
Stresses and strains are related to each other via the constitutive equations er(x) = L(x)e(x) + a(x) or 
e(x) = M(x)er(x) +/3(x), where L(x) and M(x) = L(x) _1 are the known phase stiffness and compliance 
fourth-order tensors, and the common notation for contracted products has been employed: Le = Lijki^kl- 
/3(x) and a(x) = — L(x)/3(x) are second-order tensors of local eigenstrains and eigenstresses. In particular, 
for isotropic constituents the local stiffness tensor L(x) is given in terms of the local bulk modulus fc(x) 
and the local shear modulus /x(x), and the local eigenstrain /3(x) is given in terms of the bulk component 
/3o(x) by the relations: 

L(x) = (dk,2 f j,) = dk(x)N 1 + 2n(x)N 2 , /3(x) =/3 (x)<5, (2.1) 

Ni = S ® S/d, N2 = I Ni (d = 2 or 3); 8 and I are the unit second-order and fourth-order tensors, 
and £§> denotes tensor product. For the fiber composites it is the plane-strain bulk modulus fcp] - instead 
of the 3-D bulk modulus k[ 3 ] - that plays the significant role: k[ 2 ] = k[g\ +/X[3]/3, /X[ 2 ] = (im. We introduce 
a comparison body, whose mechanical properties g c (g = L,M, a,/3,f) denoted by the upper index c 
and L c , M c will usually be taken as uniform over w, so that the corresponding boundary value problem 
is easier to solve than that for the original body. All tensors g (g = L, M, a, (3) of material properties are 
decomposed as g = g c + gi(x) = g c + gj 7 "'(x) at x e v^ m \ The upper index ( m ) indicates the components 
and the lower index i indicates the individual inclusions; = w\v, v = Uv^ = U«j, V(x) = = 

^It is known that for 2-D problems the plane-strain state is only possible for material symmetry no lower than orthotropic 
(see e.g. [14]) that will be assumed hereafter in 2-D case. 
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J2 Vi(x), and V<*>(x) and ^(x) are the indicator functions ofv^ and v il respectively, equals 1 at x 6 
and otherwise, (m = 0, fc; fc = 1, 2, . . . , TV; i = 1,2,...). 

We assume that the phases are perfectly bonded, so that the displacements and the traction com- 
ponents are continuous across the interphase boundaries, i.e. [[<xn mt ]] = and [[u]] = on the interface 
boundary T mt where n mt is the normal vector on T mt and [[(.)]] is a jump operator. The traction 
t(x) = er(x)n(x) acting on any plane with the normal n(x) through the point x can be represented in 
terms of displacements t(x) = t(n, V)u(x) + a(x)n, where Tifc(n, V) = Lijkinj((x)d / dxi. The boundary 
conditions at the interface boundary will be considered together with the mixed boundary conditions on 
r with the unit outward normal n r 

u(x) = u r (x), xer u , (2.2) 
<r(x)n r (x) = t r (x), xeT t , (2.3) 

where T u and T t are prescribed displacement and traction boundaries such that r„ur t = T, r„nr t = 0. 
u r (x) and t r (x) arc, respectively, prescribed displacement on r„ and traction on T t ; mixed boundary 
conditions, such as in the case of elastic supports are possible. Of special practical interest are the 
homogeneous boundary conditions 

u r (x) = e r x, e r = const., x e T, (2.4) 
t r (x) = cr r n r (x), er r = const., x e T, (2.5) 

where e r (x) = i [V <E) u r (x) + (V <S> u r (x)) T ] , x £ T, and e r and er r are the macroscopic strain and 
stress tensors, i.e. the given constant symmetric tensors. We will consider the interior problem when the 
body occupies the interior domain with respect to T. 

2.2 Statistical description of the composite micro structure 

It is assumed that the representative macrodomain w contains a statistically large number of realiza- 
tions a (providing validity of the standard probability technique) of inclusions V{ £ of the constituent 
7j( fe ) (i = 1,2,...; k = 1, 2, . . . , TV). A random parameter a belongs to a sample space A, over which a 
probability density p(a) is defined (see, e.g., [15, 16]). For any given a, any random function g(x, a) (e.g., 
g = V, V^ k \ cr, e) is defined explicitly as one particular member, with label a, of an ensemble realization. 
Then, the mean, or ensemble average is defined by the angle brackets enclosing the quantity g 

(g)(x)= f S {x,a)p{a)da. (2.6) 
J A 

No confusion will arise below in notation of the random quantity g(x, a) if the label a is removed. One 
treats two material length scales (see, e.g, [17]): the macroscopic scale L, characterizing the extent of w, 
and the microscopic scale a, related with the heterogeneities Vi. Moreover, one supposes that applied field 
varies on a characteristic length scale A. The limit of our interests for both the material scales and field 
one is 

L>A>a. (2.7) 

All the random quantities under discussion are described by statistically inhomogeneous random fields. For 
the alternative description of the random structure of a composite material let us introduce a conditional 
probability density (p(vi, Xj|i>i, xi, . . . , v n , x n ), which is a probability density to find the «-th inclusion 
with the center x, in the domain Vi with fixed inclusions vi, . . . ,v n with the centers xi , . . . , x n . The 
notation tp(vi, Xj|; v±, xi, . . . ,w„,x„) denotes the case x, ^ xi, . . . ,x„. We will consider a general case of 
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statistically inhomogeneous media with the homogeneous matrix (for example for so-called Functionally 
Graded Materials (FGM)), when the conditional probability density is not invariant with respect to 
translation: ip(vi,x.i + x|vi,xi, . . . ,w„,x„) ^ f{vi, Xj|i>i, xi + x, . . . , v„,x„ + x), i.e. the microstructure 
functions depend upon their absolute positions (see e.g. [18]). In particular, a random field is called 
statistically homogeneous in a narrow sense if its multi-point statistical moments of any order are shift- 
invariant functions of spatial variables. Of course, tp(vi, Xj|; i>i, Xi, . . . , u n ,x n ) = for values of x, lying 
inside the "excluded volumes" UwJJ,, (since inclusions cannot overlap, m = 1, . . . ,ri), where v^i 3 v m with 
indicator function is the "excluded volumes" of Xj with respect to v m (it is usually assumed that 
v mi = v m)> and f( v h x »l; v ii x i> • ■ ■ ) v n, *n) — > <p(vi,Xi) as |xj - x TO | ► oo, m = 1, ... ,n (since no long- 
range order is assumed). <p(vi,x) is a number density, = n^ k \x) of component 9 at the point 
x and = c^(k) is the concentration, i.e. volume fraction, of the component u, € v^ k ' at the point 
x: C W(x) = (yW)(x) = u i? 7.( fc )(x), Ui = mem (k = 1, 2, . . . , N; i = l,2,...), C ( 0) (x) = 1 - (V)(x). 
The notations ((.))(x) and ((.)|«i,Xi; . . . ;u n ,x n )(x) will be used for the average and for the conditional 
average taken for the ensemble of a statistically inhomogeneous field X = (u,) at the point x, on the 
condition that there are inclusions at the points xi,...,x n and Xj =/= Xj if i ^ j (i, j = l,...,n). 
The notations ((.)|; x^; . . . ;u„,x„)(x) arc used for the case x ^ The notation ((■)) i (x) at 

x 6 d, C v^ k ' means the average over an ensemble realization of surrounding inclusions (but not over the 
volume m of a particular inhomogeneity, in contrast to ((-))(i)) & t the fixed V{. 

2.3 General integral equation for composites of any structure 

In the framework of the traditional scheme, we introduce a homogeneous "comparison" body with 
homogeneous moduli L c , and with the inhomogeneous deterministic transformation field /3 c (x) and body 
force f c (x) (and with solution er°, e°, u° to the same boundary- value problem). For all material tensors 
g (L, M, (3, a, f) the notation gi(x) = g(x) — g c is used. 

Then, substituting the constitutive equation and the Cauchy equation into the equilibrium equa- 
tion, we obtain a differential equation with respect to the displacement u which can be reduced to a 
symmetrized integral form after integrating by parts (see, e.g, [19] and Chapter 7 in Ref. [10]) 

e(x) = e°(x) + f U(x - y)x(x)dy 

J W 

+ f VG(x-y)f 1 (y)dy+ f VG(x - s)r(s)n(s)ds, (2.8) 
Jw Jr 

where r(x) = Li(y)[e(y) — /3(y)] — L c /3 1 (y) is called the stress polarization tensor, and the surface 
integration is taken over the external surface T with the outer normal n(s)J_r of the macrodomain 
w C R d , and the integral operator kernel U is an even homogeneous a generalized function of degree — d 
defined by the second derivative of the Green tensor G: Uijki(x) = [VjV;Gifc(x)] n.-.n^y the parentheses 
in indices mean symmetrization, and G is the infinite-homogeneous-body Green's function of the Navier 
equation with an elastic modulus tensor L c defined by 



V jl/i [V <g> G(x) + (V <g> G(x)) T ] | = -<M(x), 



(2.9) 



and vanishing at infinity (|x| — > oo), <5(x) is the Dirac delta function and S is the unit second order tensor. 
The deterministic function £ (x) is the strain field which would exist in the medium with homogeneous 
properties L c and appropriate boundary conditions (see, e.g, [20]): 



S 9 ( x ) = / Gi(p, g )( x -s)ti(s) -Ui(s)L? jk[ G kMl (x-s)nj(s) ds+ / G i(Ptq) (x - y)f-(y)dy, (2.10) 



On the thermo-elastostatics of heterogeneous materials. I 



5 



which conforms with the stress field er°(x) = L c e°(x) — /3 c (x), t = t a c ■ n(s). The representation 
(2.10) is valid for both the general cases of the first and second boundary value problems as well as for 
the mixed boundary-value problem (see for references [10]). In particular, for the conditions (2.2), (2.4) 
and /3 c ,f c = the right-hand-side integral over the external surface in (2.10) can be considered as a 
continuation of £ r (x), x e T„ e T i.e. (2.4), into w as the strain field that the boundary condition 
(2.2), (2.4) would generate in the comparison medium with homogeneous moduli L c . For simplicity we 
will consider only internal points x £ w of the microinhomogencous macrodomain w at sufficient distance 
from the boundary 

a< |x-s|, VseT, (2.11) 

when the validity of Eq. (2.10) takes place except in some "boundary layer" region close to the surface 
s e T where some boundary data [u(s),t(s)] (if they are not prescribed by the boundary conditions) will 
depend on perturbations introduced by all inhomogeneities, and, therefore e°(x) = e°(x, a). 

It should be mentioned that for the constant gravitation loads ff = p c gi with a constant mass 
density p c and a constant gravitation field gi, as well as for a centrifugal load fi = gijXj with the matrix 
gij = const., the volume integral in Eq. (2.10) can be transformed into a surface integral (see e.g. [20]). 
Consider next the construction of the regularization of generalized function of the type of derivatives of 
homogeneous regular function we will use a scheme proposed by Gel'fand and Shilov [21] according to 
which the tensor U(x) is split into two parts 

U(x) =U s (x)+U / (x), (2.12) 

where U 5 (x) = <5(x)U , (U = const.) is a singular function associated with some infinitely small 

f 

exclusion region and U / (x) = l/r- d U (n), (x = ra, r = |x|) is a formal function. Both terms on the 
right-hand side of (2.12) depend on an exclusion region being prescribed, while their sum, being the 
left-hand side of (2.12), is defined uniquely (see, e.g., [10]). 



3. Random structure composites 

3.1 Known general integral equations 

To avoid much mathematical manipulations, we will consider in this subsection the case of a pure 
mechanical loading f,/3 = and L c = I/ ' =const. Then for a finite number of inhomogeneities totally 
placed in the macrodomain vo, Eq. (2.8) is reduced to simplified equation 

e(x) = £°(x)+ / U(x-y)L 1 (y)e(y)dy, (3.1) 

J w 

where e°(x) =constant for the homogeneous boundary conditions (2.4) being considered. In early mi- 
cromechanical research Eq. (3.1) was also exploited for the limiting case of a statistically homogeneous 
field of an infinite number of inhomogeneous in the whole space w = R d . This unjustified generalization 
leads to well known convergence difficulties because U(x — y) is homogeneous generalized function of de- 
gree —d and the integral in (3.1) is only conditionally convergent that gives no meaningful results without 
additional justification for the mode of integration employed. The nature of the conditional convergence 
is that the value of an integral taken over an infinite domain depends on the shape of this domain 
in specifying the limit process. Lipinski et al. [22] reduced the system (2.8) and (2.10) (at f, (3 = 0) 
to an equation analogous to (3.1); the subsequent convergence difficulty was overcome by the use of a 

self-consistent approach which is equivalent in fact to the termination of the constituent U (n) (2.12) 
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leading to the known convergence troubles. Ju and Tseng [23, 24] also used Eq. (3.1) and eliminated the 
difficulties induced by the dependence of a conditionally convergent integral on the shape of integration 
domain w by the use of an assumption of this shape (see for details Chapter 7 in [10]). Fassi-Fehri et al. 
[25] postulated the size of integration domain in Eq. (3.1). For the purely mechanical loading (f , (3 = 0) 
Eq. (2.8) formally coincides with the analogous relations obtained by Lipinski et al. (1995) by the use of 
Green's functions for a bounded domain w when the surface integral in the right-hand side (2.8) vanishes. 
Its implementation is not trivial because the finite-body Green's function G w is generally not known, and 
replacing G w by G (3.1) if w is large enough leads to well known convergence difficulties as discussed 
in detail by Willis [13]. These difficulties mentioned above and some other ones can be avoided by a few 
ways. 

One way of modifying such a conditionally convergent integral resulted by the long-range interactions 
is the so-called method of normalization (or renormalization, in analogy to its use in quantum field theory) 
achieved by subtracting from (3.1) the conditionally convergent behavior which is asymptotically closed 
to U(x — y)Li(y)e(y) at |x — y| — > oo. The renormalization procedure consists in subtraction out of the 
conditionally convergent term by making use of an expression that has the identical convergence properties 
as the term presenting the difficulty, but which also has a limiting value that is known. The renormalization 
method was systematically developed for rheological problems [26] , for conductivity of random suspensions 
[27], and for elasticity of composites with spherical particles [28] in terms of perturbations introduced 
by the dipole strengths of inhomogeneities. However, the correct choice of a rcnormalizing quantity is 
not always straightforward that initiated Willis and Acton [29] (see also [30]) to propose an alternative 
method for obtaining a convergent integral expressed through the Green function. Rigorous justification 
of the last approach was proposed by O'Brien [31] by applying the divergence theorem to the boundary 
integral in the equation analogous to Eq. (2.8) that leads to (for /3,f = 0, e° = (e)) 

e(x) = e°+ [ U(x- y)[L 1 (y) £ (y) - (L l£ )]dy, (3.2) 

J w 

where the operation of "separate" integration of slow U(x — y) and fast Li(y)e(y) variables will be 
analyzed with Eq. (3.9). Comparing of non renormalized (3.1) and renormalized (3.2) equations, McCoy 
[19, 32] suggested that one can formally remove the conditionally convergent term appearing in the former 
by simply setting them equal to zero. There arc well known "noncanonical rcgularizations" proposed by 
Kroner [33] (see also Kroner [34, 35]) and attributed to Kanaun [36] (see also [37]) 

J U(x - y)h dy = 0, or J T(x - y)h dy = L c h, (3.3) 

and 

j U(x - y)h dy = M c h, or J T(x - y)h dy = 0, (3.4) 

for the first and the second boundary-value problem, respectively; h is an arbitrary constant symmetric 
second-order tensor, and the integral operator kernel, 

T(x - y) = -L c [L5(x - y) + U(x - y)L c ] , (3.5) 

called the Green stress tensor (see [35]) is defined by the second derivative of the Green tensor G (2.9). 
According to (3.5), each of the relations is a consequence of the other from the same pair, either (3.3) 
or (3.4). In light of the note by McCoy [19, 32], the essence of the "noncanonical regularization" (3.3) 
[and (3.4)] used in Eq. (3.1) can be considered as some sort of the renormalization method. In actuality 
this regularization is an intuitive introduction of an operation of generalized functions U and T on a 
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constant symmetric tensor h = const, (e.g. h = (Lie) =const., see for comparison Gel'fand and Shilov 
[21]); Buryachenko [38] proved that the correctness of this regularization is questionable. 

The renormalizcd quantity in the O'Brien method (3.2) is obtained directly from the macroscopic 
boundary term at the homogeneous boundary conditions (2.4) that simultaneously defines both the 
advantage with respect to the classical renormalization method (see e.g. [29]) and the fundamental limi- 
tation of possible generalizations of the proposed approach to both the functional graded materials and 
inhomogeneous external loading [because (Lie) =const. in Eq. (3.2)]. 

The mentioned deficiency of Eq. (3.2) could be avoided by a centering method systematically devel- 
oped by Shermergor [12] also for statistically homogeneous media. This method was generalized in Ref. 
[39] and justified by Buryachenko [38] as applied to the FGMs. Indeed, the centering method consists in 
subtracting from Eq. (2.8) their statistical average yielding (at f, (3 = 0) 

e(x) = (e)(x) + f U(x - y)[L 1 (y)e(y) - {Ue){y)]dy. (3.6) 

The statistical averages (e)(x), (Lie)(y) T^const. contained in Eq. (3.2) allows one to apply this equation 
for the analyses of nonlocal effects appearing in both the FGMs and a case of inhomogeneous external 
loading (see [40]; and Chapter 12 in [10]). 

The original purpose of the renormalizing term was only to provide the absolute convergence of 
the integral in Eq. (3.1) that is archived by long-range behavior of the function U(x — y)(Lie) at 
|x — y| — > co. However, the same term is exploited in a short-range domain |x — y| < 3a in the vicinity 
of the point x g w. This term is defined only by the average tensor (Lie) and ignores any available 
microtopological information (possible nonellipsoidal and laminated structure of inhomogeneities, the 
shape and size of excluded volume, radial distribution function, orientation, etc.). We will demonstrate 
in the next Subsection that Eqs. (3.2) and (3.6) can be improved by the use of a new renormalizing term 
unequally determined from rigorous averaging scheme and directly dependent on the microtopological 
information mentioned above. 

3.2 New general integral equations 

As noted, the prospective centering method is based on subtracting from both sides of Eq. (2.8) their 
statistical averages. Now we will center Eq. (2.8) by the use of statistical averages presented in a general 
form (U(x — y)g)(y), i.e. from both sides of Eq. (2.8) their statistical averages are subtracted 

e(x) = (e)(x)+ f ((VG(x- y )fi))(y)dy+Xf 

+ f ((U(x- y ){Li(y)[e(y)-/3(y)]-L c /3 1 (y)}))(y)dy, (3.7) 

where one introduces a centering operation 

((U(x - y)g»(y) = U(x - y)g(y) - (U(x - y)g)(y), (3.8) 
and in the right-hand-side of Eq. (3.7), the integral over the external surface T 

«VG(x - s){Li(s)[e(s) - /3(b)] - Li/3 c }))(s)n(s)ds + e°(x, a) - (e°)(x), (3.9) 

can be dropped out, because this tensor vanishes at sufficient distance x from the boundary T (2.11). This 
means that if |x— s| is large enough for Vs 6 T, then at the portion of the smooth surface ds ps |x— s\ d ~ 1 du> s 



8 



V. A. Buryachenko 



with a small solid angle du) s the tensor VG(x — s)|x — s| d_1 depends only on the solid angle u) s variables 
and slowly varies on the portion of the surface ds; in this sense the tensor VG(x — s) is called a "slow" 
variable of the solid angle u) s while the expression in curly brackets on the right-hand-side integral of 
Eq. (3.9) is a rapidly oscillating function on ds and is called a "fast" variable. Therefore we can use 
a rigorous theory of "separate" integration of "slow" and "fast" variables, according to which (freely 
speaking) the operation of surface integration may be regarded as averaging (see for details, e.g., Ref. 
[41] and its applications Shermergor [12]). If (as we assume) there is no long-range order and the function 
(p(vj,Xj\;vi,Xi) — ip(vj,Xj) decays at infinity (as |xj — Xj| — > oo) sufficiently rapidlj{| then it leads to a 
degeneration of both the surface integral (3.9) and the summand £°(x, a) — (e°)(x). 

In order to express Eq. (3.7) in terms of stresses we use the identities: 

Li(e-/3) = -L c Mi<t, (3.10) 
e = [M c <r + /3 C ] + [Mkt + ^J. (3.11) 

Substituting (3.10) and (3.11) into the right-hand-side and the left-hand-side of (3.7), respectively, and 
contracting with the tensor L c gives the general integral equation for stresses 

<x(x) = (<t)(x) + f [«r(x - y)77))(y) + «L c VG(x - y)fi»(y)]4y+l£ (3.12) 

where we define 

II = J «L c G(x-s)L c 77))(s)n(s)ds, (3.13) 

r/(y) = M 1 (y) < r(y)+/3 1 (y). (3.14) 

If we assume no long-range order, then the tensor X T a is degenerated and can be dropped. The tensor rj is 
called the strain polarization tensor and is simply a notational convenience. In (3.14) Mi(y) and y9 1 (y) 
are the jumps of the compliance M( fe ) and of the eigenstrain (3^ inside the component t/ fe ' (k = 0, . . . , N) 
with respect to the constant tensors M c and (3 C , respectively. 

For convenience of the forthcoming presentation we will recast Eq. (3.12) in another form, for which 
we introduce the operation 

A x (x) = A(x) - (A) (x) (3.15) 

for the random function A (e.g. A = rr, e, rj, f ) with statistical average in the matrix (A)o(x). Then Eq. 
(3.12) can be rewritten in the form 

< t(x) = ( < t)(x)+ [ [((r(x-y)r, 1 ))(y) + «L c VG(x-y)fJ))(y)]dy. (3.16) 

The general integral equations (3.7), (3.12) and (3.16) are new and proposed for statistically inho- 
mogeneous media for the general case of inhomogencity of tensors /3°(x), j3^°\x), f c (x), f^(x) for both 
the general case of the first and second boundary value problems as well as for the mixed boundary- 
value problem. In some particular cases these equations are reduced to the known ones that will be 
demonstrated in Subsection 3.3. 

§ Exponential decreasing of this function was obtained by Willis [42] for spherical inclusions; Hansen and 
McDonald [43], Torquato and Lado [44] proposed a faster decreasing function for aligned fibers of circular cross- 
section. 
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All integrals in Eqs. (3.7), (3.12) and (3.16) converge absolutely for both the statistically homogeneous 
and inhomogeneous random fields X ol inhomogeneities. Indeed, even for the FGMs, the term ((r(x — 
y) r ?))(y) (3-12) is of order 0(|x - y|~ 2d+2 ) as |x - y| -> oo, and the integrals in Eqs. (3.7), (3.12) 
with the kernels U and T, respectively, converge absolutely. In a similar manner, the integrals with the 
body force density converge absolutely. In fact, the kernel VG(x — y) is of order 0(\x — y|~ d+1 ) as 
|x — y| — > oo, and the term ((VG(x — y)f 1 ))(y) tends to zero with |x — y| — > oo (x e Vi, y € Vj) as 0(|x — 
y\~ d+1 ){ip(vj , Xj\; Vi, Xj) — (p(vj,Xj)}. For no long-range order assumed, the function ip(vj ,Xj\;vi,Xi) — 
(p(vj,Xj) decays at infinity sufficiently rapidly and guarantees an absolute convergence of the integrals 
involved. Therefore, for x € w considered in Eqs. (3.7), (3.12) and (3.16) and removed far enough from 
the boundary T (a <C |x — s|, Vs e T). the right-hand side integrals in (3.7), (3.12) and (3.16) do not 
depend on the shape and size of the domain w, and they can be replaced by the integrals over the whole 
space R d . With this assumption we hereafter omit explicitly denoting R d as the integration domain in 
the equations 

e(x) = ( £ >(x) + J [«U(x - y)r))(y) + «VG(x - y)fi»(y)]dy, (3.17) 

<x(x) = (<x}(x)+! [((r(x-y)r 7 ))(y) + ((L c VG(x-y)f 1 ))(y)]dy, (3.18) 

<x(x) = ( CT )(x) + J [((r(x-y)r 7 1 ))(y) + ((L c VG(x-y)fl))(y)]dy (3.19) 

corresponding to Eqs. (3.7), (3.12), and (3.16); in a similar manner the domain w can be replaced by the 
whole space R d and omitted. Thus, there arc no difficulties connected with the asymptotic behavior of 
the generalized functions VVG and T decaying at infinity as |x — y\~ d , and there is no need to postulate 
cither the shape or the size of the integration domain w [45] or to resort to either regularization [33] , [36] 
or renormalization ([28], [32], see also Willis [13]) of integrals which are divergent at infinity. The rigorous 
mathematical analysis of correctness of these above mentioned methods is beyond the purpose of the 
current article. Nevertheless, it should be noted that the disruption of statistical homogeneity of media 
often leads to additional difficulties, whose resolution by these known mentioned approaches appears to 
be questionable (see for details Subsection 3.2.3 in Ref. [10]). 

3.3 Some particular cases 

The subsequent analysis of Eqs. (3.17)-(3.19) can be done for the comparison medium with any 
elastic modulus L c , which necessarily leads to some additional assumptions for the structure of the strain 
fields in the matrix (see for details Chapter 8 in [10]). Equations (3.17)-(3.19) are much easier to solve 
when they contain the stress-strain fields only inside the heterogeneities. There arc two fundamentally 
different approaches to ensuring it. 

In the first one we postulate 

L C = L (0) . (3.20) 

Then the integrands with the arguments y in Eqs. (3.17)-(3.19) vanish at y e v^. However, it does not 
guarantee a protection from the necessity of estimation of stress-strain distributions in the matrix in the 
general cases of both the inhomogeneous inclusions and inhomogeneous boundary conditions. Fortunately, 
this domain of the matrix is only located in the vicinity of a representative inhomogeneity v q (see for 
details [10]). 

In the second one we choose L c quite arbitrarily, and analyze Eq. (3.18) [Eq. (3.17) can be considered 
analogously]. Equation (3.18) being exact for any (rj)o(x) can be simplified with the additional assumption 
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that the strain polarization tensor in the matrix tj(x), (x G v^ ') coincides with its statistical average in 
the matrix 

77(x) = (r,) (x), xe«(°). (3.21) 

In so doing, the assumption (3.20) is more restricted in the sense that the assumption (3.20) yields the 
assumption (3.21) (the opposite is not true) and, moreover, in such a case the exact equality tj(x) = 
(tj)o(x) ee 0, x G u (0) holds. 

Equations (3.17)-(3.19) contain the general representations for statistical averages such as, e.g., 
((U(x — y)g))(y), ((r(x — y)g))(y) (g = Li£,Li/3, rj). We will consider the particular cases of these 
equations obtained for the different particular approximations of statistical averages (U(x — y)x) (y) and 
(r(x — y)rj)(y) in Eq. (3.7) and (3.12), respectively. Such an analysis will be performed for Eq. (3.17) 
[Eqs. (3.18) and (3.19) can be considered analogously] for no body forces acting and purely mechanical 
loading i.e. 

f (x) = 0, /3(x) = (3.22) 
The deterministic analog of the mentioned approximations can be presented in the following forms 

/u(x-y)p(y)Vi(y)dy = «,U(x - x i )(p) (3.23) 

/u(x-y)p(y)V;(y)dy - ^T|(x - x<)(p) (i)> (3.24) 

where p(y) (y £ t)j) is some deterministic function, Vi is some representative fixed heterogeneity, and the 
tensors 

T £ Cx-x )-i -( v i)~ lp i forxGUi, . . 

l^x xj-| ( _ rl/u(x _ y) ^ (y)dy forx$? ^, l^&J 

have analytical representations for ellipsoidal inclusions in an isotropic matrix (see for reference [10]), 
and Pj = -/ U(x - y)^(y)dy ee SiM (0) =const. (for Vx G v t ) is defined by the Eshelby [46] tensor S. ; . 
Obviously that the equalities (3.23) and (3.24) are only asymptotically fulfilled at |x — Xj| — ► oo, and it 
is possible to propose a formal counterexample where an error of both approximations (3.23) and (3.24) 
equals infinity, e.g. if p(x) ^ and (p)(j) = (x G u,). The most popular approximation (3.23) (which 
is simultaneously the most crude) was implicitley used by many authors (see for early references, e.g., 
[47]) including implied exploiting of Eq. (3.6) for obtaining of some sort of the centered Eq. (3.2) (see 
[30]). A quantative analysis of results obtained by the use of the representations (3.23) and (3.24) will be 
performed in an accompanied paper by Buryachenko [48]). 

Substitution of the random analog [e.g., when p(x) is replaced by g(x, a) = r(x, a), tj(x, a)] of the 
approximation (3.23) into Eq. (3.17) for statistically homogeneous media subjected to the homogeneous 
boundary conditions (2.4) yields the known Eq. (3.2) with the renormalizing term obtained by O'Brian [31] 
at /3(x) = through the homogeneous boundary conditions. Moreover, for the mentioned homogeneous 
boundary conditions and statistically homogeneous media (n^(x) = =const., q = 1,...,N), the 
approximations (3.23) and (3.24) leads to an identical result reducing Eq. (3.17) to (3.2). This statement 
holds if we prove that contributions made to Eq. (3.17) by the renormalizing terms (3.23) and (3.24) are 
identical for any macrodomain x G w (/3(x) ee 0): 

N 

/ U(x-y)dy(L ie ) / T e g (x - Xq )v q n^ (x q )(L l£ ) q dx q (3.26) 

J W i J W 
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For justification of the equality (3.26), it should be mentioned that for uniform distribution of inclusion 
centers x g , all volume of the domain w in the right-hand side of Eq. (3.26) is uniformly covered by the 
moving ellipsoids v q . Then any point in the domain x G w in the right-hand side integral (3.26) is covered 
by the same number k of the ellipsoids v q with homogeneous strain polarization tensor i~(y) = (Lie) 
(y £ v q): an d, therefore, the integral over the covered domain w on the right-hand side of Eq. (3.26) is 
equal (within some probability factor) to k integrals over domain w in the left-hand side. Therefore, in 
the case r(y) =const. inside moving inhomogcneity y <G v q , both approximation (3.23) and (3.24) reduce 
Eq. (3.17) to the known one (3.2). However, a condition of homogeneity r(y) =const. at y £ v q is fulfilled 
only for homogeneous ellipsoidal inhomogencities in the framework of an additional hypothesis of effective 
field homogeneity according to which each inclusion is located inside a homogeneous socalled effective 
field (see also [48]). An abandonment from effective field hypothesis leads with necessity to inhomogcneity 
of the stress-strain fields inside the inhomogencities that can tend to the different predictions of effective 
moduli based on Eqs. (3.2) and (3.17) even for both the statistically homogeneous media and homogeneous 
boundary conditions (see for details [48]). This difference is a result of insensitivity of the renormalizing 
term U(x — y)(Lie) [obtained at the approximation (3.23)] in the correct equation (3.2) to the details 
of heterogeneities of the stress-strain fields inside the inclusions, while a corresponding term (U(x — 
y)Lie)(y) [which is exact and obtained without approximations neither (3.23) nor (3.24)] of Eq. (3.17) 
explicitly depends on the mentioned field inhomogcneity. 

It should be mentioned that the equality used in Eq. (3.26) (e.g., g = Lie) 

N 

(g> = 5> (9) <g>, (3.27) 

q=l 

is only fulfilled for statistically homogeneous media subjected to the homogeneous boundary conditions; 
here the summation in the right-hand side is performed over the volume of the representative inclusions 
v q G (q = 1, . . . , N). If any of these conditions is broken then it is necessary to consider two sorts 
of conditional averages (see for details [10]). At first, the conditional statistical average in the inclusion 
phase (g)^(x) = (gV) l ' q \x.) (at the condition that the point x is located in the inclusion phase x G v^) 
can be found as (gV)^ q \x) = (V^ (x)) _1 (gV^ qS) )(x). Usually, it is simpler to estimate the conditional 
averages of these tensors in the concrete point x of the fixed inclusion x G v q : (g|u g ,x g )(x) = (g) g (x). 
Although in a general case 

N N 

(g)(x) = ]T C ^(g) (<?) (x) £ £ c W(gK,x g >(x), (3.28) 

where v q G v^ q ' , it can be easy to establish a straightforward relation between these averages for the 
ellipsoidal inclusions v q with the semi-axes a g = (a q , . . . , a q ) T . Indeed, at first we built some auxiliary 
set Wq(x) with the boundary <9t^(x) formed by the centers of translated ellipsoids v q (0) around the fixed 
point x. We construct v q (x) as a limit v® — > v q (x) if a fixed ellipsoid Vk is shrinking to the point x. Then 
we can get a relation between the mentioned averages [x = (xi, . . . , a;<i) T ]: 

(g) (9) (x)=/ rc (,) (y)(gK,y)(x) dy. (3.29) 

J«J(x,a,) 

Formula (3.29) is valid for any material inhomogcneity of inclusions of any concentration in the macrodomain 
w of any shape (if w^(x) C w). Obviously, the general Eq. (3.29) is reduced to Eq. (3.27) for both the 
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statistically homogeneous media subjected to homogeneous boundary conditions and statistically homo- 
geneous fields g (e.g., g = er, e). 

Thus, we have performed a qualitative competitive analysis of both the known (3.2), (3.6) and new 
(3.17), (3.19) general integral equations. Quantitative correlation of some estimations obtained at the 
bases of these equations is presented by Buryachenko [48]. Interested readers are referred to the book 
by Buryachenko [10] for detailed comparison of Eqs. (3.2) and (3.6) with the related equations and 
approaches. It should be mentioned, that the particular cases of Eq. (3.6) were also widely used (either 
explicitly or implicitly) by other authors (see, e.g., [15], [16], [49-52]). However, Eq. (3.6) were obtained in 
the mentioned papers by a centering method based on subtracting from Eq. (3.1) [rather than from Eq. 
(2.8)] their statistical average obtained in the framework of implicit use of the asymptotic approximation 
(3.23) (although this approximation was not indicated). 

Lastly, we will consider the field X bounded in one direction such as a laminated structure of some 
real FGM (see [1], [2]). Then the surface integral (2.8) over a "cylindrical" surface (with the surface area 
proportional to p = |x — s|) tends to zero with |x — s| — > oo as p~ d+2 simply because the generalized 
function VG(x — s) is an even homogeneous function of order — d + 1. Therefore, for infinite media the 
surface integral (2.8) vanishes, and Eq. (2.8) can be rewritten as 



e(x) = e°(x) + / VG(x - y)f 1 (y)ciy + / U(x- y){L 1 (y)[e(y) - /3(y)] - L^^dy, (3.30) 



In so doing the integrals from body forces in Eqs. (3.30) and (3.31) only conditionally converge for a 
general case of a bounded function fi(y); because of this, for the function f(y) we will assume decay 
at infinity |x — y| — > oo no less then 0(|x — y| _/3 ), (/? > 0) guaranteeing the absolute convergence 
of body force integrals in Eqs. (3.30) and (3.31). Clearly in the considered case of X bounded in one 
direction, Eqs. (3.30) and (3.31) are exact, and the right-hand-side integrals in (3.30) and (3.31) converge 
absolutely. Eq. (3.30) was used by Torquato [17], [53] for the particular case (3.30) with homogeneous 
boundary conditions (2.4) and for the inclusion field X with a constant concentration of inclusions within 
an ellipsoidal domain included in the infinite matrix. Although Eqs. (3.6), (3.17) and (3.19) are more 
complicated then Eqs. (3.30) and (3.31), nevertheless they provide practical advantages because their 
integrands decay at infinity faster then the integrands involved in Eq. (3.30) and (3.31). 

It should be noted that Eqs. (3.17)-(3.19) exploiting the infinite-homogeneous-body Green functions 
were obtained from Eqs. (3.7), (3.12), and (3.16), respectively, at sufficient distance x from the boundary 
T (2.11), and, therefore, they can not be used for analysis of boundary layer effects (e.g. free edge effect). 
For this class of problems, the Green functions for finite domains seem more prospective (see, e.g., [15], 
[54] and Chapter 14 in [10]). In such a case, the appropriate general integral equations generalizing Eqs. 
(3.17)-(3.19) to the finite domains can be obtained in a straightforward manner (see for details Chapter 
14 in [10]). However, more detailed considerations of boundary layer effects are beyond the scope of the 
current study and will be analyzed in other publications. 

4. Random structure composites with long-range order 

Localized Eqs. (3.17) and (3.18) were obtained in the framework of no long-range order assumption 
when the integrand in curly brackets decays at infinity |x — s| — > oo sufficiently rapidly. Now we relax 



or, alternatively, in terms of stresses 




(3.31) 
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this assumption and for the sake of definiteness, we will consider some conditional averages of the surface 
integral in Eq. (3.18) 

(/SKxi; • • ■ ; v«,xn)(x) = J {(r r (x - s)T7 X (s)|; «i,x i; . . . ; v n ,x n )(s) - (r r (x - s)j 7 1 (s))(s)}ds, (4.1) 

where x £ V\, . . . , w„, (n = 1,2, . . .), x T and T r (x — s) = — L c VG(x — s)L c . The asymptotic behavior 
of the integrand in curly brackets in Eq. (4.1) as |x— s| — *■ oo can be estimated by the use of representation 
of the solution (j7 1 |;w 1 ,x 1 ; . . . ;u„,x„)(s) by the successive approximation method (see for details [10]). 
Then 

(T r (x - a)* 1 (s) | ; v 1: x i; . . . ; v n , x„) (s) - (T r (x - s)^ 1 (s)) (s) 

-* (r^X - s)?7 1 (s)) i (s)[(p(u ! ;,X ? ;|wi,Xi, . . . ,W„,X„) - ^(Wj,Xj)] 
n 

+ 0(r- 2d+1 )^(r ? 1 ),(x J )^(w 4 ,x> 1 ,x 1 ,...,t.„,x„), (4.2) 
i=i 

where s £ Vi, r = min \xj — s|, (j = 1, . . . , ri) and the terms in Eq. (4.2) of order 0(r~ 3d+1 ) and higher 
order terms are dropped. The contribution of terms in (4.2) proportional to 0{r~ 2d+1 ) in the integral 
(4.1) are degenerated at |xj — s| — > oo and Eq. (4.1) can be simplified 

(£^|«i,xi; . . . ;w„,x„)(x) = J (T r (x - s)r7 1 (s)) l (s)[^(w i , Xi|ui, xi, . . . ,u„,x„) - ^(wj, Xj)]efe. (4.3) 

If boundary conditions are applied for which (er)(s) and, therefore, (i7 1 )j(s) vary linearly (or higher) 
with s and [ip(vi,Xi\vi,Xi, . . . , u„,x„) — y)(vj,Xi)] does not decay sufficiently rapidly as — s| — * cxd 
(j = 1, . . . , ft ) (long-range order) then the integral (4.3) may be divergent. We will consider the interesting 
practical case of a random structure composite described as either a triply or a doubly periodic in the 
broad sense random field X. 

Namely, it is now assumed that the representative macrodomain w contains a statistically large 
number of realizations of ellipsoidal inclusions u,; £ C R d (i = 1,2,.. .) with identical shape, orienta- 
tion and mechanical properties. The composite material is constructed using the building blocks or cells: 
w = Uf2 m , v m C f2 m . We consider a composite medium with each random realization of particle centers 
distributed at the nodes of some spatial lattice A with the nodes corresponding to j = 1, . . . , n parti- 
cles in each cell i7 m . Suppose ej (i = 1, . . . , d) are linearly-independent vectors, so that we can represent 
any node m £ A of both the triply and doubly periodic lattice as 



d-l 



£[£' 



i + f d (m d )e d \, (4.4) 

3=1 *=1 3=1 »=1 

where m-' = (m\, . . . , m 3 d ) are integer-valued coordinates of the node m-? in the basis ej which are equal 
in modulus to |e»| , and f 3 d {mi) — f 3 d {mi + 1) ^ const., (i = 1, . . . ,d). In the plane /(m<j) = const, the 
composite is reinforced by periodic arrays A md of inclusions in the direction of the ei axis and the e^-i 
axis. The type of the lattice A md is defined by the law governing the variation in the coefficients m, (i = 
1, 2), and also by the magnitude and orientation of the vectors e, (i = 1, d— 1). In the functionally graded 
direction the inclusion spacing between adjacent arrays may vary {f 3 d (nid) — fd(md + 1) ^ const.). 
For a doubly-periodic array of inclusions in a finite ply containing 2m' + 1 layers of inclusions we have 
f 3 (rrid) = at \md\ > m l ; in more general case of doubly periodic structures f 3 (rrid) ^ at m<j — > ±oo. 
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To make the exposition more clear we will assume that the basis e t is an orthogonal one and the axes 
(i = 1, . . . , d) are directed along axes of the global Cartesian coordinate system (these assumptions are 
not obligatory). 

Let V x be a "moving averaging" cell with the center x and characteristic size ay = \/v, and let 
for the sake of definitencss £ be a random vector uniformly distributed on V x whose value at z E V x is 
(pg(z) = 1/V X and </?g(z) = otherwise. Then we can define the average of the function g with respect 
to translations of the vector £ for each random realization of the function g 

(g)x(x-y) = =- J g(z-y)dz, xGlV (4.5) 

Among other things, "moving averaging" cell V x can be obtained by translation of a cell f2i and can vary 
in size and shape during the motion from point to point. Clearly, contracting the cell V x to the point x 
occurs in passing to the limit (g) x (x — y) — > g(x — y). To make the exposition more clear we will assume 
that V x results from £li by translation of the vector x — xp; it can be seen, however, that this assumption 
is not mandatory. 

The surface integral (4.3) can be eliminated in the equation related with (3.12) by "centrification" 
achieved by subtracting from both sides of Eq. (3.12) their averages over the moving averaging cell V x 
(4.5). In so doing the average operator (4.5) introduced for a deterministic function g(y) should be recast 
for random function g(y) by the use of a previous estimation of a statistical average (g)(z — y): 



(g) x (x-y) = =j- / (g)(z-y)dz, xefi, (4.6) 
Then Eq. (3.12) is reduced to 

«r(x) = <<x)(x) + f [«r(x-y)r,)) x (y) + <(L c VG(x - y)f 1 )) x (y)] dy + ({X r J x , (4.7) 

J w 

where 

((z£)) x = J ((T r (x - s)»7 1 (s))) x (s)[(/3(f l ,x l |ui,xi, . . .,v n , x„) - (p(vi,Xi)]ds. (4.8) 

is a centered surface integral and one introduces a new centering operation (((-))) x such as, e.g., 

«T(x - y)i?» x (y) = T(x - y)T,(y) - (T(x - y)n) x (y). (4.9) 

For the analysis of integral convergence in Eq. (4.7) at |x — y| — » oo, we assume that »7 1 (y) can be 
regarded as a constant, equal to the value at Xj = y.j, and may thus be taken outside the averaging 
operation (({■)}} x . Then we expand T(z — y) (z € Vx) in a Taylor series about x and integrate term by 
term over the cell V x with the center x 

T(z - y) = T(x - y) + (z - x)Vr(x - y) + i(z - x) <g> (z - x)VVr(x - y) . . . , (4.10) 

(D x (x-y) = r(x-y) + -^-/ (z-x)®(z-x)rfzVVr(x-y).... (4.11) 
The similar expansion can be performed for the tensor VG(x-y) that leads to 

«r)) x (x-y) = (z-x)®(z-x)dzVVr(x-y)+.... (4.12) 

^ VxJ]/ x 

<(VG)) x (x-y) = --i-/ (z-x)®(z-x)rfzVVVG(x-y)+.... (4.13) 
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As is evident from Eq. (4.12), the tensor ((r)) x (x— y) is of order 0(ay |x— y| d 2 ) with the dropped terms 
in Eq. (4.12) being of order 0(a^|x — y|~ d ~ 4 ) and higher order terms. Then the absolute convergence of 
volume integral (4.7) is assured because at sufficient distance x from the boundary r as |x — y| — > oo the 
integration over y can be carried out independently for both ((U)) x (x — y) (the function of the "slow" 
variable x — y) and Tj(y) (the function of "fast" variable y), and therefore the volume integral converges 
absolutely. In a similar manner the term ((r r (x — s)rj 1 (s))) x (s) in the surface integral (4.8) is of order 
0(ay |x— y| _d_1 ), and the surface integral vanishes at |x — s| — > oo , s G V if (<x)(s) and, therefore (r] 1 )j(s) 
grows with s slower then 0(|x — s| 2 ~^) (}3 = const. > 0). For the same reason the volume integral with 
the integrand ((L c VG(x — y)fi)) x (y) converges absolutely for bounded functions fi(y). 

By this means the locality principle exists in the new Eq. (4.7) for the case of long-range order 
composites being considered if the average stress (<r)(s) grows with s slower then 0(|x — s\ 2 ~@) (j3 = 
const. > 0). 

5. Conclusion 

As noted in Introduction, the final goals of micromechanical research of composites involved in a 
prediction of both the overall effective properties and statistical moments of stress-strain fields are based 
on the approximate solution of the exact initial integral equation (2.8). Absolute convergence of integrals 
in Eq. (2.8) is provided by different versions of the centering procedure performed for the different cases of 
the boundary conditions, microtopology of composite material, and accompanied assumptions. It leads to 
the different general Eqs. (3.2), (3.6), and (3.17) (and their stress-state analogs) called the initial integral 
equations which have the different renormalizing terms. Then, considering some conditional ensemble 
averages of the general equations either (3.2), (3.6), or (3.17) yield the infinite hierarchy of equations. 
These truncated hierarchies of equations are solved as a system of coupled equations. One starts with 
the last members of these hierarchies, the ones which have the most inclusions held fixed, because these 
equations does not depend on the others. The fields so obtained give the previous terms in the next 
equations up the hierarchies. One continues step by step up the hierarchies until the unconditionally 
averaged fields arc finally obtained. However, these standard procedures (differing by both the numbers of 
coupled equations and assumptions exploited for their solutions) have fundamentally diverse backgrounds 
defined by the features of the renormalizing terms in Eqs. (3.2), (3.6), and (3.17). 

These features of the initial integral equations are fundamental for subsequent solving the truncated 
hierarchy involving a rearrangement of each appropriate equation before it is solved. The most successful 
rearrangement are those which make the right-hand side of the coupled equations reflect the detailed 
corrections to that basic physics. So, the advantage of Eq. (3.2) with respect to Eq. (3.30) (even for the case 
when Eq. (3.30) is correct) is explained by faster convergence of corresponding integrals as |x — y| — > oo. 
The centering method realized at the obtaining of Eq. (3.2) subtracts the difficult state at infinity from 
the equation, i.e. roughly speaking the constant forcc-dipolc density expressed through an alternative 
technique of the Green's function. This dictates the fundamental limitation of possible generalization of 
Eq. (3.2) to both the FGMs and inhomogeneous boundary conditions. The mentioned deficiency of Eq. 
(3.2) was resolved by Eq. (3.6) which renormalizing term provides an absolute convergence of the integral 
in Eq. (3.6) at |x — y| — ► oo for general cases of the FGMs. However, the same term in Eq. (3.6) is used 
in a short-range domain |x — y| < 3a in the vicinity of the point x G w. A fundamental deficiency of 
Eq. (3.6) is a dependence of the renormalizing term U(x — y)(Lie)(y) [obtained in the framework of the 
asymptotic approximation (3.23)] only on the statistical average (Lie)(y) rather than (U(x — y)Lie)(y) 
in Eq. (3.17). What seems to be only a formal trick [abandoning the use of the approximations of the 
kind (3.23) and (3.24)] is in reality a new background of micromcchanics (3.17)-(3.19) [which does not 
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use an approximation of the kind either (3.23), (3.24), or (3.3) as in Eqs. (3.2) and (3.6)] yielding the 
discovery of fundamentally new effects even in the theory of statistically homogeneous media subjected 
to homogeneous boundary conditions (see for details the accompanied paper by Buryachenko [48]). 
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